In this lab, we will explore some of R’s functionality with spatial data. For some of the longer examples, it is highly recommended to use R’s scripting functions to run the examples to save on re-typing.

Before starting the lab, you will need to set up a new folder for your working directory. Go to your geog6000 folder now and create a new folder for today’s class called lab08. The following files will be used in this lab, all available on Canvas:

Climate dataset for Western North America: WNAclimate.csv
Temperature dataset for Oregon in shapefile format: oregon.zip
New York state polygon data in a shapefile: NY_Data.zip
Digital elevation map from Switzerland: swiss_dem.grd
A NetCDF file of global monthly air temperature: air.mon.ltm.nc

You will need to download these files from Canvas, and move them from your Downloads folder to the datafiles folder that you made previously. Make sure to unzip the zip files so that R can access the content. Note that on Windows, you will need to right-click on the file and select ‘Extract files’.

You will also need to install the following packages:

sp
maptools
rgdal
raster
RColorBrewer
classInt
ggmap

With all the examples given, it is important to not just type in the code, but to try changing the parameters and re-running the functions multiple times to get an idea of their effect. Help on the parameters can be obtained by typing help(functionname) or ?functionname.

Spatial Objects

`SpatialPoint`* Objects

Read in the WNA Climate dataset from the csv file WNAclimate.csv, and use the names() function to see the columns in the dataframe:

wnaclim <- read.csv("../datafiles/WNAclimate.csv")
names(wnaclim)

## [1] "WNASEQ"   "LONDD"    "LATDD"    "ELEV"    
## [5] "totsnopt" "annp"     "Jan_Tmp"  "Jul_Tmp"

We now create a simple SpatialPoint object of the location of the sites in the dataframe using the SpatialPoints() function. Note that we first create a coordinate reference system (CRS) for the data (here, simple longitude/latitude on a WGS84 sphere).

library(rgdal)
llCRS <- CRS("+proj=longlat +ellps=WGS84")

## Warning in showSRID(uprojargs, format = "PROJ",
## multiline = "NO"): Discarded datum Unknown based
## on WGS84 ellipsoid in CRS definition

wnaclim.sp <- SpatialPoints(cbind(wnaclim$LONDD, wnaclim$LATDD), 
                            proj4string = llCRS)
summary(wnaclim.sp)

## Object of class SpatialPoints
## Coordinates:
##                 min  max
## coords.x1 -138.0000 -100
## coords.x2   29.8333   60
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=WGS84 +no_defs]
## Number of points: 2012

There are some standard functions to obtain basic data and metadata about any Spatial* object:

bbox(wnaclim.sp)
proj4string(wnaclim.sp)

## Warning in proj4string(wnaclim.sp): CRS object has
## comment, which is lost in output

coordinates(wnaclim.sp)

The original dataframe also contains various other parameters (elevation, climate, site name) in columns that we can associated with the SpatialPoints, by creating a SpatialPointsDataFrame.

wnaclim.spdf <- SpatialPointsDataFrame(cbind(wnaclim$LONDD, 
                                             wnaclim$LATDD), 
                                       wnaclim, proj4string = llCRS)
summary(wnaclim.spdf)

Now we have a spatial object with associated values, we can make some plots. Using R’s basic plot() function with any spatial data will simply plot the spatial geometry:

plot(wnaclim.spdf)

The sp package has a function to make thematic maps. For this we need to specify the Spatial* object and the column name or attribute to use. Note that both this and the plot() function will adapt to different spatial data types.

spplot(wnaclim.spdf, "Jan_Tmp", colorkey=TRUE)
spplot(wnaclim.spdf, "Jul_Tmp", colorkey=TRUE)

As the attributes for these locations are held as a data frame, we can use all the standard R functions with them. For example, to calculate the mean and standard deviation of January temperatures:

mean(wnaclim.spdf$Jan_Tmp)

## [1] -3.375149

sd(wnaclim.spdf$Jan_Tmp)

## [1] 7.382896

And to plot a histogram of annual precipitation:

hist(wnaclim.spdf$annp, main="W.NAm Annual Precipitation", breaks=20)

`SpatialPolygon`* Objects

BLAKE - THIS SECTION CAN BE REMOVED Polygons can be created using a combination of the Polygon() function and the SpatialPolygon() function. Here, we create a simple, rectangular polygon and plot it as an overlay over western North America.

The polygon is built using the Polygon() function — note that the first and last coordinates must be the same to close the polygon. The next function
The Polygons() function builds a list of polygons — this allows a set of non-contiguous polygons to be collected together
The SpatialPolygons() function then converts this to a Spatial* object, which allows us to use this with other spatial data

Sr1 = Polygon(cbind(c(-130,-110,-110,-130,-130), c(30,30,50,50,30)))
Srs1 = Polygons(list(Sr1), "s1")
wusa <- SpatialPolygons(list(Srs1), proj4string=llCRS)
plot(wnaclim.spdf)
plot(wusa, add=T, col=3, lwd=2)

We can now use this to select only points within the polygon. This uses the same syntax as R’s indexing, where here it select all rows or locations that fall within the polygon:

wnaclim.sub = wnaclim.spdf[wusa,]
plot(wnaclim.spdf)
plot(wusa, add=TRUE)
plot(wnaclim.sub, pch=16, col=2, add=TRUE)

BUT NOT THIS BIT However, by far the easiest way to get polygon data into R is to export it from a GIS package. Shapefiles are one of the most commonly used formats, and the rgdal package includes a function (readOGR()) for reading these directly into R. The zip file NY_data.zip contains shape files for part of the state of New York. Unzip this in your current working directory, then read in the data as follows:

NY8 <- readOGR("../datafiles/NY_data/NY8_utm18.shp")

## Warning in OGRSpatialRef(dsn, layer, morphFromESRI
## = morphFromESRI, dumpSRS = dumpSRS, : Discarded
## datum D_unknown in CRS definition: +proj=utm
## +zone=18 +ellps=WGS84 +units=m +no_defs

class(NY8)
summary(NY8)

This creates a SpatialPolygonDataFrame — a polygon object with spatial attributes and an associated data frame. As the shapefile had an associated file with projection information (a *.prj file), the projection is automatically set in R.

proj4string(NY8)

## Warning in proj4string(NY8): CRS object has
## comment, which is lost in output

We can now plot the shapefile. As before, we use the attributes as a data frame. If, for example, we want to calculate the total population:

sum(NY8$POP8)

## [1] 1057673

And, we can again use spplot() to make a thematic map by specifying specify which column or attribute we want to use:

spplot(NY8, "POP8")

If we are only interested in a subset of the data, we can use the attributes to select polygons. The data frame contains a column with the district associated with each polygon. If we only want the subset corresponding to Syracuse, we can extract this as follows:

Syracuse <- NY8[NY8$AREANAME == "Syracuse city", ]
plot(Syracuse)
spplot(Syracuse, "POP8")

Or we can select by an attribute. For example to extract only the polygons with over 15% population ages 65 and over (this is the variable ‘PCTAGE65P’):

Pop65p <- Syracuse[Syracuse$PCTAGE65P > 0.15, ]
plot(Pop65p)
spplot(Pop65p, "PCTAGE65P")

The full set of vector formats that can be read by readOGR() can be seen by typing ogrDrivers() at the console.

`SpatialGrid`* Objects

As with the other types of objects, SpatialGrid objects exist with and without associated data frames. The first type (SpatialGridDataFrames) are used for storing raster data (remote sensed images, gridded climate fields, etc). The second type are less useful as they contain no data, only the grid topology.

The maptools package includes a function (readAsciiGrid()) for reading raster files in the ESRI ArcInfo format. The file swiss_dem.grd contains a raster data set, with a digital elevation model of Switzerland, and can be read in as follows:

library(maptools)

## Checking rgeos availability: TRUE

swiss.dem <- readAsciiGrid("../datafiles/swiss_dem.grd", colname='elev')
spplot(swiss.dem, "elev")

As this raster images has a number of missing values (outside Switzerland), we can convert it to a SpatialPixelDataFrame, where the grid topology is retained, but pixels with missing values are deleted.

## Convert to SpatialPixelDataFrame
swiss.dem2 <- as(swiss.dem,'SpatialPixelsDataFrame')
summary(swiss.dem2)

Note that a lot of R’s function for working with gridded data has been moved to the raster package - examples of the use of this are given below.

Writing out spatial data

The writeOGR() allows to write spatial data back out as a shapefile. To write out the SpatialPoint* object we made earlier:

writeOGR(wnaclim.spdf, dsn="wnaclim.shp", 
         layer="wnaclim", driver="ESRI Shapefile")

Map projections

R has functions for reprojecting spatial data using the PROJ4 libraries, and functions from the rgdal package. Try the following code to reproject a set of sites in Oregon, and the state polygon boundary.

ortann <- readOGR("../datafiles/oregon/oregontann.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/u0784726/Dropbox/Data/devtools/geog6000/datafiles/oregon/oregontann.shp", layer: "oregontann"
## with 92 features
## It has 4 fields

orotl <- readOGR("../datafiles/oregon/orotl.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/u0784726/Dropbox/Data/devtools/geog6000/datafiles/oregon/orotl.shp", layer: "orotl"
## with 36 features
## It has 1 fields

proj4string(ortann) <- CRS("+proj=longlat")
proj4string(orotl) <- CRS("+proj=longlat")

In this code, we read in the data from shape files (points and polygons), then assign a geographical projection (longitude/latitude) to both objects, using the CRS() function. There are several ways to do this. The above code simply sets the projection system to longitude/latitude. An alternative and easier way is to use EPSG codes. These are a standard set of geographic and Cartesian projections, each with a numeric code. The codes can be found here. The code for a standard long/lat projection, with a WGS84 ellipsoid and datum is 4326, so we could replace the last two line in the previous chunk of code with:

proj4string(ortann) <- CRS("+init=EPSG:4326")
proj4string(orotl) <- CRS("+init=EPSG:4326")

To reproject, we first create a projection object that contains the parameters for the new projection as name/value pairs For example ‘+proj=aea’ means set the map projection to Albers Equal Area. Other parameters include:

The name of the projection (+proj)
The two standard parallels (+lat_1 +lat_2)
The coordinates of the projection center (+lon_0 +lat_0)
The ellipsoid (+GRS80)
The units of distance (+units)

aea.proj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 
  +lon_0=-110 +ellps=GRS80 +datum=NAD83 +units=m"

Finally, we use the spTransform() function to do the reprojection. This requires a Spatial* object, and our projection system parameters, which are converted to a coordinate reference system using CRS().

orotl.proj <- spTransform(orotl, CRS(aea.proj))
ortann.proj <- spTransform(ortann, CRS(aea.proj))

If you now compare the coordinates of the projected and unprojected objects, you will see the difference:

summary(ortann)
summary(ortann.proj)
par(mfrow=c(1,2))
plot(orotl)
plot(ortann, col=2, pch=16, add=TRUE)
plot(orotl.proj)
plot(ortann.proj, col=2, pch=16, add=TRUE)

There are many other projections (and parameters) available in the PROJ4 library, and more details and examples can be found here.

Visualizing spatial data

The basic plotting function provided with the sp package, spplot() provides a simple way to map out data with a simple color scale (e.g. spplot(NY8[,"POP8"])). This has several arguments that can improve on the basic maps. The following examples demonstrate the use of some of these, and require two extra packages to be installed: classInt (provides methods for selecting data intervals for colors) and RColorBrewer (provides better color palettes).

SpatialPointDataFrame Objects

This next example uses a shapefile of country borders downloaded from Natural Earth, who provide free global vector and raster data for map making. We will start by loading this:

countries <- readOGR("../datafiles/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/u0784726/Dropbox/Data/devtools/geog6000/datafiles/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp", layer: "ne_50m_admin_0_countries"
## with 241 features
## It has 94 fields
## Integer64 fields read as strings:  POP_EST NE_ID

plot(countries)

Now we’ll plot the Western North America data set. The additional arguments are:

main: adds a title
edge.col: is the color used for the border of the points
sp.layout: includes other layers to add to the plot. The order they are included here will define the order in which they are added to the map
colorkey: when TRUE this forces a continuous color scale.

spplot(wnaclim.spdf,"Jul_Tmp", main="W.NAm July Temp. (deg C)", 
       edge.col=1, sp.layout=list(countries), colorkey=TRUE)

We can now change the palette used for the color scale, by choosing one of the RColorBrewer palettes. There are several of these, split in sequential, diverging and categorical types. These were developed by Cynthia Brewer specifically for map visualization and include palettes designed for color-blindness, projection and printing. Full details can be at the ColorBrewer website. We’ll start by loading the package and showing all available palettes:

library(RColorBrewer)
display.brewer.all()

The names of the palettes are given on left. We will used a diverging palette for the temperature data (feel free to try any of the others as well). To use this with spplot() requires two steps. First, we create the new palette, specifying the number of intervals and the name of the palette required. (Note that we reverse this to get red = warm and blue = cold.) Then we add this to the spplot() function using the col.regions argument. We also need to specify the number of cuts - which is \(k-1\), where \(k\) is the number of intervals in the palette.

my.pal <- rev(brewer.pal(n = 9, name = "RdYlBu"))
spplot(wnaclim.spdf,"Jul_Tmp", main="W.NAm July Temp. (deg C)",
       sp.layout=list(countries), colorkey=TRUE, edge.col=1, 
       col.regions=my.pal, cuts=8)

By default, RColorBrewer will only use discrete color palettes. You can, however, create a continuous palette using a combination of brewer.pal() and colorRampPalette() as follows. This second function takes the ColorBrewer palette and interpolates it to 100 colors:

col.levs=100
my.pal.cont = colorRampPalette(my.pal)(col.levs)
spplot(wnaclim.spdf,"Jul_Tmp", main="W.NAm July Temp. (deg C)",
       sp.layout=list(countries), colorkey=TRUE, edge.col=1, 
       col.regions=my.pal.cont, cuts=col.levs-1)

SpatialPolygonDataFrame Objects

We’ll now replot the Syracuse shapefile with a better color palette. We’ll add an additional function classIntervals(), which can be used to change the intervals in a color palette. By default, these are even splits (so each interval covers the same range). If your variable is skewed, then other splits may be better. The style argument allows to choose different types of split. Here, we’ll use percentiles, which acts in a similar way to a histogram normalization of the variable. Note that we now have three lines of code to make the plot:

Make the color palette
Create the intervales
Make the plot. Note the use of the at argument, which sets the breaks to those defined by classIntervals()

library(classInt)
my.pal <- brewer.pal(n = 7, name = "OrRd")
breaks.qt <- classIntervals(Syracuse$POP8, n = 6, style = "quantile")
spplot(Syracuse, "POP8", col.regions = my.pal, col = "transparent", 
       at = breaks.qt$brks)

SpatialGridDataFrame Objects

And finally an example with a gridded dataset. This is the Swiss DEM used above. We’ll also load a shapefile with the Switzerland borders:

swiss.bord = readOGR("../datafiles/borders/borders.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/u0784726/Dropbox/Data/devtools/geog6000/datafiles/borders/borders.shp", layer: "borders"
## with 4 features
## It has 1 fields

my.pal <- rev(brewer.pal(n = 9, name = "Greens"))
breaks.qt <- classIntervals(swiss.dem$elev, n = 8, style = "quantile")

## Warning in classIntervals(swiss.dem$elev, n =
## 8, style = "quantile"): var has missing values,
## omitted in finding classes

spplot(swiss.dem, "elev", main="Switzerland Elevation",
       sp.layout=list(swiss.bord), col.regions = my.pal, at = breaks.qt$brks)

The raster package

The raster package offers a array of functions for dealing with gridded data, including the ability to read from many widely used file formats, including remote sensing images (e.g. GeoTiffs) and NetCDF and HDF formats. We will use it here to work with gridded air temperature data (air.mon.ltm.nc) from the [NCAR NCEP reanalysis project][ncarID]. This is the long term means for each month from 1981-2010. The file has 12 layers (one per month) and one variable (air).

Raster objects

Start by installing and loading the package, the use the raster() function to read from the file.

library(raster)
r = raster("../datafiles/air.mon.ltm.nc")

Note that we have only read the first layer (January). R will tell you that it loaded the variable called air. To avoid this message you can specify this directly, which is important for files containing multiple variables:

r = raster("../datafiles/air.mon.ltm.nc", varname="air")
r

## class      : RasterLayer 
## band       : 1  (of  12  bands)
## dimensions : 73, 144, 10512  (nrow, ncol, ncell)
## resolution : 2.5, 2.5  (x, y)
## extent     : -1.25, 358.75, -91.25, 91.25  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : /Users/u0784726/Dropbox/Data/devtools/geog6000/datafiles/air.mon.ltm.nc 
## names      : Monthly.Long.Term.Mean.Air.Temperature.at.sigma.level.0.995 
## z-value    : 0000-12-30 
## zvar       : air

r is a raster object and has information about grid spacing, coordinates etc. Note that the description of the object tells you whether it is held in memory (small raster files) or on disk. We can make a simple plot:

plot(r, main="NCEP NCAR January LTM Tair")

And you should be able to see the outline of the continents. By default, NCAR NetCDF files have longitudes running from 0 to 360. We can convert this to the more commonly used (and UK-centric :) -180 to 180 by the function rotate(). We can also use a different color palette and add the country shapefile we loaded earlier:

r = rotate(r)
my.pal <- brewer.pal(n = 9, name = "OrRd")
plot(r, main="NCEP NCAR January LTM Tair", col=my.pal)
plot(countries, add=TRUE)

The function cellStats() can be used to calculate most summary statistics for a raster layer. So to get the mean global temperature (and standard deviation):

cellStats(r, mean)
cellStats(r, sd)

Values can be extracted from individual locations (or sets of locations) using extract(). This can take a set of coordinates in matrix form, or use a Spatial* object. To get the January temperature of Salt Lake City:

extract(r, cbind(-111.9,40.76))

By default this gives you the value of the cell in which the point falls. The value can equally be estimated by bilinear interpolation from the four closest cells with method='bilinear':

extract(r, cbind(-111.9,40.76), method='bilinear')

We created a SpatialPoints object earlier with the location of samples in Western North America (wnaclim.spdf). We can now use this, and the raster layer to get the January temperature for all locations:

extract(r, wnaclim.spdf, method='bilinear')

If we want to use only a subset of the original raster layer, the function crop() will extract only the cells in a given region. This can be defined using another raster object or Spatial* object, or more simply by defining an extent object:

myext = extent(c(-130,-60,24,50))
r = crop(r, myext)
plot(r, main="NCEP NCAR January LTM Tair", col=my.pal)
plot(countries, add=TRUE)

Raster Stacks

A useful extension to the basic raster functions is the use of stacks. These are a stack of raster layers which represent different variables, but have the same spatial extent and resolution. We can then read in and store all 12 months from the NetCDF file, and then work with this. We read these in with stack() and crop them as before:

myext = extent(c(-130,-60,25,50))
r.stk = rotate(stack("../datafiles/air.mon.ltm.nc", varname="air"))
r.stk = crop(r.stk, myext)
r.stk

By typing the name of the stack object, we can see that this has 12 layers, each with 280 cells and the extent, etc. The names attributed to each layer are often unreadable, so we can add our own names:

names(r.stk) <- paste("TAS",month.abb)
names(r.stk)

Now any of the functions we used previously with one layer will be used across all layers. The plot() function returns a grid with one plot per layer. Setting the zlim argument ensures that all figures use the same range of colors:

plot(r.stk, col=my.pal, zlim=c(-35,35))

Adding a shapefile (or other spatial information) is a little more complex. We create a simple function to plot the country borders, then include this as an addfun in the call to plot():

addBorder = function(){ plot(countries, add=TRUE) }
plot(r.stk, zlim=c(-35,35), col=my.pal, addfun=addBorder)

The cellStats() function now returns the mean (or other statistic) for all layers, allowing a quick look at the seasonal cycle of average air temperature.

tavg = cellStats(r.stk, mean)
plot(1:12, tavg, type='l', xlab="Month", ylab="Avg T (C)")

And we can do the same for an individual location using extract():

slc.tavg = extract(r.stk, cbind(-111.9,40.76), method='bilinear')
plot(1:12, slc.tavg, type='l', xlab="Month", ylab="Avg T (C)")

This same approch allows you to extract pixels by polygon overlays. Let’s re-read the global dataset, and try extracting values for a given country. First, we’ll extract China’s polygon as a new SpatialPolygon, then use this in the extract() function to get all the pixels within the borders

r = rotate(raster("../datafiles/air.mon.ltm.nc"))

## Warning in .varName(nc, varname, warn = warn): varname used is: air
## If that is not correct, you can set it to one of: air, valid_yr_count

china.sp = subset(countries, NAME == "China")
extract(r, china.sp)

As this function can be used with a set of polygons, the output is in a list. We can use extract the first value of this list, and make a quick histogram:

china.tjan = extract(r, china.sp)[[1]]
hist(china.tjan)

The extract() function also takes an argument fun. This allows you to calculate a summary statistic for each set of pixels that is extracted (i.e. one per polygon). Here, we’ll use this with the countries SpatialPolygon to get an average value of January temperature. We add this back as a new column in the countries object, and then plot it:

countries$Jan_Tmp = extract(r, countries, fun = mean)[,1]
spplot(countries, "Jan_Tmp", main = "Country average January temperature")

Other spatial functionality

R and KML files

R can write out KML files for use with Google Earth. A few things to note here:

The function is kmlPoints() from the maptools package (there are also functions for KML polygons and lines)
The function requires a Spatial* object (here a SpatialPointsDataFrame), and the name of an output file
Optionally we can provide a layer name (kmlname) and a vector of names for each point (name)
The optional parameter description associates some descriptive text with each point, including HTML tags.
- Here we create a series of descriptions that include the January and July temperatures, together with text labels.
- Note the use of the paste() function to glue together different text strings (which can include numeric vectors)

The second function writes out the points and descriptions into a KML file in the current working directory. Once run, find this file and open it in Google Earth. You can replace the icon using any of the icons found here.

description = paste("<b>TJan:</b>", wnaclim.spdf$Jan_Tmp, 
                    "<br><b>TJul:</b>", wnaclim.spdf$Jul_Tmp)
kmlPoints(wnaclim.spdf,"test.kml", kmlname="WNAclim",
          name=wnaclim$WNASEQ, description=description, 
          icon="http://www.google.com/mapfiles/kml/paddle/wht-diamond.png")

Here’s a second example, with polygon data (the Syracuse polygons). Note that as these are in a UTM projection, we need to reproject to latitude/longitude coordinates, then can write out to a KML file:

ll.proj <- "+proj=longlat +ellps=WGS84"
Syracuse.proj <- spTransform(Syracuse, CRS(ll.proj))

kmlPolygons(Syracuse.proj, kmlfile = "syr.kml", 
            col="white", border="black", lwd=3)

Appendix: The tmap package

There are a couple of very good packages that extend the basic visualization of spplot(). Here, we’ll take a quick look at tmap for producing thematic maps. Another useful, but more complex package is leaflet, which allows you to produce interactive maps that can be embedded in html pages. To run these examples, you will need to have loaded the data in the examples earlier in the lab, as well as having installed tmap.

tmap works by building a series of layers and map geometry and elements. We start by using tm_shape() to identify the spatial object to be used, and then geometries are added, including filled polygons, borders, legends, etc.

library(tmap)

Example 1: Syracuse data

We’ll start by making some maps with the Syracuse dataset we made earlier. First, let’s make a simple map showing the polygon outlines using tm_borders():

tm_shape(Syracuse) + tm_borders()

## Warning in sp::proj4string(obj): CRS object has
## comment, which is lost in output

The function tm_fill() will then fill these using one of the variables in the Syracuse data set (POP8). Note that this automatically adds a legend within the frame of the figure:

tm_shape(Syracuse) + tm_borders() + tm_fill("POP8")

## Warning in sp::proj4string(obj): CRS object has
## comment, which is lost in output

The color scale can be changed by setting the palette argument in tm_fill(). This includes the ColorBrewer scales described above,and the different intervals. For example, to use the ‘Greens’ palette with percentile breaks:

tm_shape(Syracuse) + tm_borders() + 
  tm_fill("POP8", palette = "Greens", style = "quantile")

## Warning in sp::proj4string(obj): CRS object has
## comment, which is lost in output

Other map elements can be added. Here we add a longitude/latitude graticule with tm_graticules(), a north arrow and a line of text with the date the map was made.

tm_shape(Syracuse) + tm_graticules(col = "lightgray") + tm_borders() + 
  tm_fill("POP8", palette = "Greens", style = "quantile") + 
  tm_compass(position = c("left", "bottom")) + 
  tm_credits("2019-10-19", position = c("right", "top"))

## Warning in sp::proj4string(obj): CRS object has
## comment, which is lost in output

Interactive

tmap also has a function to produce a simple, interative map. This is done by changing the mode of operation from plot to view. The following code set the mode to interactive, makes a simple plot of the Syracuse data, and resets the mode to static plots. Note that these plots will be interactive in RStudio, and can also be embedded in the html files produced in R Markdown.

tmap_mode("view")
names(Syracuse)
tm_shape(Syracuse) + tm_borders() + tm_fill("Cases", palette = "Greens")
tmap_mode("plot")

Example 2

We’ll next make some maps with the Western North American site data. Individual symbols can be plotted on a color scale using tm_symbols.

tm_shape(wnaclim.spdf) + tm_symbols(col="Jan_Tmp")

## Warning in sp::proj4string(obj): CRS object has
## comment, which is lost in output

This takes the same arguments as tm_fill() for the color palette. We’ll use a red to blue color scale from RColorBrewer. The minus sign before the palette name reverses the order of the colors. As there is a large amount of overlap between the sites, we also add an alpha level to make the symbols transparent.

tm_shape(wnaclim.spdf) + 
  tm_symbols(col="Jan_Tmp", alpha = 0.5, palette = "-RdBu")

## Warning in sp::proj4string(obj): CRS object has
## comment, which is lost in output

We’ll next add country boundaries from the Natural Earth shapefile loaded earlier. Note that as this is a different spatial object, we have to use tm_shape() a second time to reference this, then use tm_borders() to add the lines.

tm_shape(wnaclim.spdf) + 
  tm_symbols(col="Jan_Tmp", alpha = 0.75, palette = "-RdBu") + 
  tm_shape(countries) + tm_borders(col="gray")

## Warning in sp::proj4string(obj): CRS object has
## comment, which is lost in output

We can also use this

tm_shape(wnaclim.spdf) + 
  tm_symbols(col="Jan_Tmp", alpha = 0.75, palette = "-RdBu") + 
  tm_shape(countries) + tm_borders(col="gray") + tm_style("cobalt")

## Warning in sp::proj4string(obj): CRS object has
## comment, which is lost in output

Example 3

In the final example, we’ll make figures using the global air temperature dataset. Start by re-reading the data (we also change the name of the layer, to reduce the amount of typing later on…).

r = raster("../datafiles/air.mon.ltm.nc", varname="air")
r = rotate(r)
names(r) <- "jan_tmp"
proj4string(r) <- CRS("+init=epsg:4326")

Now let’s make a plot using tm_raster(), which again takes similar options for color palettes. We’ll also add the country borders.

tm_shape(r) + 
  tm_raster("jan_tmp", style="fisher", palette="-RdBu") +
  tm_shape(countries) + tm_borders()

## Warning in sf::st_is_longlat(shp2): bounding box
## has potentially an invalid value range for longlat
## data

## Warning in sf::st_is_longlat(shp): bounding box
## has potentially an invalid value range for longlat
## data

We can improve this a little by moving the color legend outside of the plotting area. We’ll increase the number of color classes to 9, and add a histogram showing the frequency of different values.

tm_shape(r) + 
  tm_raster("jan_tmp", style="fisher", palette="-RdBu", legend.hist = TRUE, n = 9) +
  tm_shape(countries) + tm_borders() + tm_layout(legend.outside = TRUE, 
                                                 legend.outside.position = "left")

## Warning in sf::st_is_longlat(shp2): bounding box
## has potentially an invalid value range for longlat
## data

## Warning in sf::st_is_longlat(shp): bounding box
## has potentially an invalid value range for longlat
## data

Many more examples of plots are available through the tmap page on CRAN here.

GEOG 6000 Lab 08 Spatial Data in R

Simon Brewer

September 01, 2020

Spatial Objects

`SpatialPoint`* Objects

`SpatialPolygon`* Objects

`SpatialGrid`* Objects

Writing out spatial data

Map projections

Visualizing spatial data

SpatialPointDataFrame Objects

SpatialPolygonDataFrame Objects

SpatialGridDataFrame Objects

The raster package

Raster objects

Raster Stacks

Other spatial functionality

R and KML files

Appendix: The tmap package

Example 1: Syracuse data

Interactive

Example 2

Example 3

GEOG 6000 Lab 08 Spatial Data in R

Simon Brewer

September 01, 2020

Spatial Objects

SpatialPoint* Objects

SpatialPolygon* Objects

SpatialGrid* Objects

Writing out spatial data

Map projections

Visualizing spatial data

SpatialPointDataFrame Objects

SpatialPolygonDataFrame Objects

SpatialGridDataFrame Objects

The raster package

Raster objects

Raster Stacks

Other spatial functionality

R and KML files

Appendix: The tmap package

Example 1: Syracuse data

Interactive

Example 2

Example 3

`SpatialPoint`* Objects

`SpatialPolygon`* Objects

`SpatialGrid`* Objects